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ON SOLVING L-SRl TRUST-REGION SUBPROBLEMS 


JOHANNES BREST, JENNIFER B. ERWAY, AND ROUMMEL F. MARCIA 


Abstract. In this article, we consider solvers for large-scale trust-region sub¬ 
problems when the quadratic model is defined by a limited-memory symmetric 
rank-one (L-SRI) quasi-Newton matrix. We propose a solver that exploits the 
compact representation of L-SRI matrices. Our approach makes use of both 
an orthonormal basis for the eigenspace of the L-SRI matrix and the Sherman- 
Morrison-Woodbury formula to compute global solutions to trust-region sub¬ 
problems. To compute the optimal Lagrange multiplier for the trust-region con¬ 
straint, we use Newton’s method with a judicious initial guess that does not 
require safeguarding. A crucial property of this solver is that it is able to com¬ 
pute high-accuracy solutions even in the so-called hard case. Additionally, the 
optimal solution is determined directly by formula, not iteratively. Numerical 
experiments demonstrate the effectiveness of this solver. 


1. Introduction 


In this article, we describe a method for minimizing a qnadratic fnnction dehned 
by a limited-memory symmetric rank-one (L-SRl) matrix subject to a two-norm 
constraint, i.e., for a given Xk, 


minimize 

pea?" 


Q{p) = g'^p + ^p'^ Bp 


subject to ||p ||2 < S, 


( 1 ) 


where g = V/(xfc), B is an L-SRl approximation to V‘^f{xk), and 5 is a given 
positive constant. In large-scale optimization, solving ([^ represents the bulk of the 
computational effort in trust-region methods. In this article, we propose a solver 
that is able to solve ([^ to high accuracy. 

High-accuracy L-SRl subproblem solvers are of interest in large-scale optimiza¬ 
tion for two reasons: (1) In previous works, it has been shown that more accurate 
subproblem solvers can require fewer overall trust-region iterations, and thus, fewer 
overall function and gradient evaluations [ZlIHli; and (2) it has been shown that 
under certain conditions SRI matrices converge to the true Hessian-a property that 
has not been proven for other quasi-Newton updates [3]. While these convergence 
results have been proven for SRI matrices, we are not aware of similar results for 
L-SRl matrices. However, we hope that this paper will facilitate the study of L-SRl 
quasi-Newton trust-region methods. 

Solving large trust-region subproblems dehned by indehnite matrices are espe¬ 
cially challenging, with optimal solutions lying on the boundary of the trust-region. 
Since L-SRl matrices are not guaranteed to be positive dehnite, additional care 
must be taken to handle indehniteness and the so-called hard case (see, e.g., mm)- 
To our knowledge, there are only two solvers designed to solve the quasi-Newton 
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subproblems to high accuracy for large-scale optimization. Specihcally, the MSS 
method [I] is an adaptation of the More-Sorensen method [TS] to the limited- 
memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton setting. More 
recently, in [1], Burdakov et ah solve a trust-region subproblem where the trust 
region is dehned using shape-changing norms. It should be noted that while the 
focus of [1] is solving trust-region subproblems dehned by shape-changing norms 
instead of the usual Euclidean two-norm, Burdakov et ah also present a trust- 
region method that is able to solve L-BFGS quasi-Newton subproblems to high 
accuracy dehned by the usual Euclidean two-norm. In this article, we present a 
method that extends what is presented in [T] to the indehnite case by handling 
three additional non-trivial cases: (1) the singular case, (2) the so-called hard case, 
and (3) the general indehnite case. We know of no high-accuracy solvers designed 
specihcally for L-SRl trust-region subproblems for large-scale optimization of the 
form ([^ that are able to handle these cases associated with SRI matrices. It should 
be noted that large-scale solvers exist for the general trust-region subproblem that 
are not designed to exploit any specihc structure of B. Examples of these include 
LSTRS [201121] and SSM (12] [Tl]. 

Methods that solve the trust-region subproblem to high accuracy are often based 
on the optimality conditions for a global solution to the trust-region subproblem 
(see, e.g.. Gay m, More and Sorensen [T8| or Gonn, Gould and Toint 0 ) , given 
in the following theorem: 

Theorem 1. Let 6 be a positive constant. A vector p* is a global solution of the 
trust-region subproblem & if and only ^/||p *||2 < d and there exists a unigue a* > 0 
such that B -\- a*I is positive semidefinite and 

{Ba*I)p* =—g and (T*(5 — ||p*|| 2 ) = 0. (2) 

Moreover, if B + a*I is positive definite, then the global minimizer is unigue. 

The MorGSorensen method [18] seeks a solution pair of the form {p*, a*) that sat- 
ishes both equations in (|^ by alternating between updating p* and a*] specihcally, 
the method hxes a*, solving for p* using the hrst equation and then hxes p*, solving 
for a* using the second equation. In order to solve for p* in the hrst equation, the 
MorGSorensen method uses the Gholesky factorization of B al] for this reason, 
this method is prohibitively expensive for general large-scale optimization when 
B does not have a structure that can be exploited. However, the MorGSorensen 
method is arguably the best direct method for solving the trust-region subproblem. 
While the MorGSorensen direct method uses a safeguarded Newton method to hnd 
a*, the method proposed in this article makes use of Newton method’s together with 
a judicious initial guess so that safeguarding is not needed to obtain a*. Moreover, 
unlike the MorGSorensen method, the proposed method computes p* by formula, 
and in this sense, is an iteration-free method. 

This article is organized in Eve sections. In the second section, we review the 
L-SRl quasi-Newton matrices how to hnd its eigenvalues. In the third section, we 
describe the proposed OBS method. Numerical results are presented in the fourth 
section, and concluding remarks are found in the hfth section of the article. 

Notation. Unless explicitly indicated, || • || denotes the vector two-norm or its 
subordinate matrix norm. The identity matrix is denoted by /, and its dimension 
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depends on the context. Finally, we assnme that all L-SRl npdates are compnted 
so that the L-SRl matrix is well defined. 


2. L-SRl MATRICES 


In this section, we review L-SRl matrices, their compact formnlation, and how 
to compnte their eigenvalues. 


Given a continuously differentiable function f{x) G 3^"" and iterates {x^}, the 
SRI quasi-Newton method generates a sequence of matrices {Bk} from a sequence 
of update pairs {{sk,yk)} where 

Sk = Xk+I - Xk and yk = Vf{xk+i) - V/(xa:), 


and V/ denotes the gradient of /. Given an initial matrix Bq, Bk is defined as 

(2/fc BkSk){yk BkSk) 


Bk+1 = Bk + 


iVk BkSk^ Sk 


(3) 


provided {yk — BkSkY'Sk 7 ^ 0. In practice, Bq is often taken to be a nonzero constant 
multiple of the identity matrix. Limited-memory SRI (L-SRl) methods store and 
use only the M most-recently computed pairs {{sk,yk)}, where M n. Often 
M may be very small (for example, Byrd et ah |1] suggest M G [3, 7]). For more 
background on the SRI update formula, please see, e.g., [I2lll5l[igiig|22ll23]. 


Compact representation. To compute a compact representation of an SRI ma¬ 
trix, we make use of the following matrices: 


Sk = [so S 2 ■■■ Sk] e gfj-x(^+b, 

Yk ^ [yo yi y2 ■■■ yk] ^ 

Furthermore, we make use of the following decomposition of S'^Yk G s)fj(^+i)x(^+i); 

Sk Yk — Lk T Dk T Rk} 


where Lk is strictly lower triangular, Dk is diagonal, and Rk is strictly upper trian¬ 
gular. We assume all updates are well-defined, i.e., {yk — BkSuYSk 7 ^ 0; otherwise, 
the update is skipped (see [m Sec. 6.2]). 

The compact representation of SRI matrices is given by Byrd et ah [H Theorem 
5.1], who showed that Bk+i in (|^ can be written in the form 

Bk+i = Bo + ^kMk^l, (4) 

where G ^ gfj{fc+i)x(fc-i-i)^ ^ diagonal matrix (i.e., Bq = 7 /, 

7 G 3?). In particular, and Mk are given by 

Tfe = Yk-BoSk and Mk = {Dk + Lk + BoSk)-\ 


Since k < M, the matrix Mk is small and can be inverted in practice. We assume 
that updates are performed so that always has full column rank and that 7 7 ^ 0 
so that Bq is invertible. 


Eigenvalues. We now review from how to compute the eigenvalues of quasi- 
Newton matrices that admit a compact representation (BIIDI). For simplicity, we 
drop the subscript k, {k < M <C n), and consider the problem of computing the 
eigenvalues of 


B = Bo + 


( 5 ) 
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where Bq = 7/, 7 G 3?. Suppose ^ = QR is the “thin” QR factorization of \I/, 
where Q G and R G is invertible. Then, 

R = 7J + QRMR^Q'^. (6) 


Let UAU'^ be the spectral decomposition of RMR^ G where U G 

gfj(fc+i)x(A:+i) jg orthogonal and A = diag(Ai,..., A^+i), with Ai < • • ■ < A^+i. We 
note that since both M and R are invertible, Aj 7 ^ 0 for i = 1,, k+1. Substituting 
this into ([^, we obtain 

B = jl + QUAU^Q^. 

Since Q and U have orthonormal columns, then T|| = QU G 3 fj«x(fc+i) g^jg^ p^g 
thonormal columns. Let P± be a matrix whose columns form an orthonormal basis 
for the orthogonal complement of the column space of T||. Then, P = [ T|| P^] g 

gfjnxn jg gach that P^P = PP^ = I. Thus, the spectral decomposition of B is given 


by 


B = PAP^ 


where A = 


Ai 0 ■ 


A -|- 7/ 0 

_0 A2 


0 7/ 


( 7 ) 


where A = diag(Ai,..., A„) = diag(Ai+7,..., Afc+1+7,7,... ,7), Ai G gfj{fc+i)x(fc+i)^ 
and A2 G gfj(fi-fc-i)x(n-fc-i)_ remaining eigenvalues, found on the diagonal of 
A2, are equal to Afc+2 = 7- (For further details, see [HllQ].) For the duration of this 
article, we assume the hrst k + 1 eigenvalues in A are ordered in increasing values, 
i.e., Ai < A2 < ... < Afc+i. Finally, throughout this article we denote the leftmost 
eigenvalue of B by Amin, which is computed as Amin = iBin{Ai,7}. 


3 . Proposed method 


The method proposed in this paper, called the “Orthonormal Basis L-SRl” (OBS) 
method, is able to solve the trust-region subproblem to high accuracy even when 
the L-SRl matrix is indehnite. The method makes use of two separate techniques. 
One technique uses (1) a Newton method to hnd a* that is initialized so its iterates 
converge monotonically to a* without any safeguarding when global solutions lie on 
the boundary of the trust region, and (2) the compact formulation of SRI matrices 
together with the strategy found in [ 2 ] to compute p* directly by formula. The other 
technique is newly proposed in this article. This technique computes an optimal 
pair (p*, a*) using an orthonormal basis for the eigenspace of B. The idea of using 
an orthonormal basis to represent p* is not new; this approach is found in [1]. In this 
manuscript, we apply this approach to the cases when B is singular and indehnite. 


We begin by providing an overview of the OBS method. To solve the trust-region 
subproblem, we hrst attempt to compute an unconstrained minimizer Pu to Q. 
If the solution exists (i.e., B is not singular) and lies inside the trust region, the 
optimal solution for the trust-region subproblem is given by p* = pu and a* = 0. 
This computation is simplihed by hrst hnding the eigenvalues of B (see 0 ) ; the 
solution Pu to the unconstrained problem is found using a strategy found in [2], 
adapted for L-SRl matrices. If ||pu|| > 5 or is not well-dehned, a global solution 
of the trust-region subproblem must lie on the boundary, i.e., it is a root of the 
following function: 


0(a) = 


1 


( 8 ) 


lb(t^)ll 

When a global solution is on the boundary, we consider three cases separately: 
(i) B is positive dehnite and ||ptj|| > 6, (ii) B is positive semidehnite, and (iii) 
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B is indefinite. We note that the so-called hard case can only occur in the third 
case. Details for each case is provided below; however, we begin by considering the 
unconstrained case. 


Computing the unconstrained minimizer. The OBS method begins by com¬ 
puting the eigenvalues of B as in Section If i? is positive dehnite, the OBS 
method computes ||pu|| using properties of orthogonal matrices. If ||pu|| < 6, then 
{p*,a*) = {pu,0). We begin by presenting the computation of ||pu||, which is only 
performed when B is positive dehnite. We include a in the derivation for complete¬ 
ness even though a = 0 when hnding the unconstrained minimizer. 

The unconstrained minimizer p^ is the solution to the hrst optimality condition 
in (|^; however, the unconstrained minimizer can also be found by rewriting the 
optimality condition using the spectral decomposition of B. Specihcally, suppose 
B = PAP^ is the spectral decomposition of B, then 

—g = {B + <jl)p = (PAP^ + al)p = P{A + al)v, 


where v = P^p. Since P is orthogonal, ||n|| = ||p||, and thus, the hrst optimality 
condition expressed in Q can be written as 

{A + aI)v = -P'^g. (9) 


Note that spectral decomposition of B transforms the hrst system in ([^ into a 
solve with a diagonal matrix in ([^. If we express the right hand side as 


P^g = [ ^1 




A 

^11 

[pisi 


g^. 


then 


ibwir = ii^(^)ir 


V (f/ii)? ], \\9±r 

'^(A* + a)2J (7 + a)2' 


( 10 ) 


Thus, the length of the unconstrained minimizer pu = p(0) is computed as ||pu|| = 
||r;(0)||, where g\i = P[g = {QUfg = {^fR-^Ufg and \\g±f = \\gf - H^ylp. 
Notice that determining ||p^i|| does not require forming p^ explicitly. Moreover, we 
are able to compute without having to compute g± = Pj^g, which requires 
computing P±, whose columns form a basis orthogonal to P\\. 

If WpuW < 6, then p* = pu and a* = 0. To compute pu, we use the Sherman- 
Morrison-Woodbury formula for the inverse of i? as in [2] , adapted from the BFGS 
setting into the SRI setting in ([^: 


p* = -^ [I - g, 


( 11 ) 


where r* = 7 . Notice that this formula calls for the inversion of -|- T^T); 

however, the size of this matrix small {{k -t-1) x (A; -|- 1) where k < M), making the 
computation practical. 

On the other hand, if ||p„|| > S, then the solution p* must lie on the boundary. 
We now consider the three cases as mentioned above. 


Case (i): B is positive definite and ||p„|| > 5. Since the unconstrained min¬ 
imizer lies outside the trust region and ||pu|| = ||p(0)||, then 0 (ct) given by (|^ is 
such that 0(0) < 0. In this case, the OBS method uses Newton’s method to hnd 
a*. (Details on Newton’s method are provided in Section 3.1 ) Finally, setting 
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T* = 7 + CT*, the global solution of the trust-region subproblem, p*, is computed 


using ( 11 ). 


Case (ii): B is singular and positive semidefinite. Since 7 7 ^ 0 and B is 

positive semidefinite, the leftmost eigenvalue is Ai = 0. Let r be the multiplicity of 
the zero eigenvalue; that is, Ai = A 2 = ... = = 0 < A^+i- For a > 0, the matrix 

(A -I- al) is invertible, and thus, ||p(o')|| in ( |Io| is well-defined for a G (0,cxd). If 
limo—).o+ 0(o') < 0, the OBS method uses Newton’s method to hnd a*. (Details on 
Newton’s method are provided in Section 3.1[ ) Setting t* = 7-(-cr*, p* is computed 
using ( 11 ). 


We now consider the remaining case: (j){a) > 0. By [HI Lemma 7.3.1], 

is strictly increasing on the interval (0, cxd). Thus, (j) can only have a root in 
the interval [0, 00 ] at cr = 0. We now show that {p*, a*) is a global solntion of the 
trust-region snbproblem with a* = 0 and 


p 


= -Btg = _p(A + a-/)tF7, 


( 12 ) 


where f denotes the Moore-Penrose psendoinverse. The second optimality condition 
holds in ([^ since a* = 0. It can be shown that the first optimality condition holds 
by using the fact that g must be perpendicular to the eigenspace corresponding to 
the 0 eigenvalne of B, i.e., {P^g)i = 0 for i = 1 ,..., r (see |IH])- 

In this subcase, the trust-region subproblem solution p* can be computed as 
follows: 


p* = -P{A + a*iyP^g 

-P\\{Ai + a*I)^P(g 
-P\l{Ai + a*I)-^P[g 
-'^R-W{Ai + a*I)^g\\ - 
-'^R-^U{Ai + a*I)-^g\\ 


J + a" 


:P±Pl9 


if a* ^ - 7 , 
otherwise 


7 -h cr" 


if a* ^-7 
otherwise. 


(13) 


which makes use of the chain of following chain of eqnalities: P±Pj^g = (/ — 
P\lP^)g = (/ — ^R~^R~^'^^)g. The actnal compntation of p* in (13) reqnires only 
matrix-vector prodncts; no additional large matrices need to be formed to find a 
global solntion of the trust-region snbproblem. 

Case (iii): B is indefinite. Since B is indehnite, Amin = iFiin{Ai, 7 } < 0- Let r 
be the algebraic mnltiplicity of the leftmost eigenvalue. For a > —Amin, (A -|- al) 
is invertible, and thus, ||p(cr)|| in ( 10 ) is well defined in the interval (—Amin, 00 ). 

If lini^_j._;^+ ^((j) < 0, then there exists a* G (—Amin, 00 ) with 0(a*) = 0 that 


can obtained as in Case (i) using Newton’s method (see Sec. 3.1). The solntion p* 
is compnted via ( 11 ) with t* = 7 -|- a*. 

If lim .+ 0(cr) > 0, then g must be orthogonal to the eigenspace associated 

^ “^min 

with the leftmost eigenvalue of B [TH]. In other words, if Amin = Ai, then {g\\)i = 0 
for z = 1,..., r; otherwise, if Amin = 7 , then || 5 fj_|| = 0. We now consider the cases 
of eqnality and inequality separately. 

If lim = 0, then a* = —Amin > 0, and a global solution of the 

■^min 

trust-region snbproblem is given by 


p* = -{B + a*iyg = -P{A + a*iyp'^ 


9- 
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As in Case (ii), p* is obtained from (13) and can be shown to satisfy the optimality 
conditions (|^. 

Finally, if lim^_ 5 .__;,^+ 0(cr) > 0, then 


lim ||p((t)|| = lim \\ — {B + a*I) 


This corresponds to the so-called hard case. The optimal solution is given by 

p* = p* + z*, where p* = — {B + cr*/)^ g, z* = aWmin, (14) 

and where Wmin is an eigenvector associated with Amin and a is computed so that 


||p*|| = 6 [IE]. As in Case (ii), we avoid forming P± using (13) to compute p*. 
The computation of Wmin depends on whether Amin is found in Ai or A 2 in ([^. If 
Amin = Ai then the hrst column of P is a leftmost eigenvector of P, and thus, Wmin 
is set to the hrst column of -F||. On other hand, if Amin = 7 , then any vector in 
the column space of P± will be an eigenvector of B corresponding to Amin- Since 
Range(P||)-*- = Range(Px), the projection matrix (/ —T||T|j^) maps onto the column 
space of Pj_. For simplicity, we map one canonical basis vector at a time (starting 
with Cl) into the space spanned by the columns of P± until we obtain a nonzero 
vector. Since dim(F||) = k + 1 ^ n, this process is practical and will result with 
a vector that lies in Range(Pj_); that is, Umin = (4 — P\\Pi^)ej for at least one j in 
{1 < j < fc -|- 2} with lluminll 7 ^ 0. (We note that both Ai and 7 cannot both be 
Amin since Ai = Ai -|- 7 and Ai 7 ^ 0 (see Section^ 


The following theorem provides details for computing optimal trust-region sub¬ 
problem solutions characterized by Theorem 1 for the case when B is indehnite. 


Theorem 2. Consider the trust-region subproblem given by 

minimize Q{p) = g^p-p^Bp subject to ||p ||2 < <5, 

2 

where B is indefinite. Suppose B = PAP^ is the spectral decomposition of B, and 
without loss of generality, assume A = diag{Xi,..., \n) is such that Amin = Ai < 
A 2 < ... < An- Further, suppose g is orthogonal to the eigenspace associated with 
•^min; Lc., g^Psj = 0 for j = l,...,r, where r > 1 is the algebraic multiplicity 
of Xmin- Then, if the optimal solution of the subproblem is with a* = —Amin? then 
the global solutions to the trust-region subproblem are given by p* = p* + z* where 
p* = — [B -\- cr*/)^ g and z* = iaWmin; where -Umin is a unit vector in the eigenspace 
associated with Amin and a = ■ Moreover, 


Q{p* ±az*) = 


-<rS. 


( 15 ) 


Proof. By [18], a global solution of trust-region subproblem is given by p* = p* + z* 
where p* = — {B -\- cr*)^ g, z* = dUmin, and a is such that ||p*|| = 5. It remains to 


show that both roots of the quadratic equation 




= 5^ are given by 


a = ±-^ 0 ^ — ||p*|p and that (15) holds. 

To see this, we begin by showing that (p*)'^ z* = 0. Let r > 1 be the algebraic 
multiplicity of Amin- Then, p* = —{B + (y*lfig = —P(A -|- a*lfiP'^g = —Pv{a*), 
where n(cr*) = (A -|- a*lfiP^g. Notice that by definition of the pseudoinverse, 
v{a*)i = 0 for i = 0 ,..., r. Since Umm is in the eigenspace associated with Amin, 
then it can be written as a linear combination of the first r columns of P, i.e.. 
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^min e:= UiPci for some {wj} G 3fJ where e* denotes the ith canonical basis 
vector. Then, 


{p*)^ z = a {p*)'^ Urain = a{Pv{a*))'^ ( ^ UiPci I = an(cr*)'^ ^ UiCi = 0, 


. ^=1 


2 = 1 


since the hrst r entries of n(cr*) are zero. Since p* is orthogonal to z*, then 


a 


n* II 2 


To see (15), consider the following: 


Q{P* ±au^in) = {P* ± aMmin)Y + ^{P* ± a U^nin)^ B {P* ± « Mmin) 

= {p* ±aUram)'^{g “ ^5' “ iP* ± OMmin)) 

= hp* ±aUram)^g -]:(T*\\p* ±aUrain\\‘^ 


T * *r 2 

= ^9P -^oS, 


( 16 ) 


Since ut 


^g = (Ei=i UiPci)^ g = (Ei=i UicfP'^) g = 0 since g is orthogonal to the 


eigenspace associated with Amin. D 


□ 


The OBS method is summarized in Algorithm 1 . 


Compute T = QR, the “thin” QR factorization (or, compute the Cholesky 
factor R of 

Compute RMRF = UAU'^ (the spectral decomposition) with 


0 ); 


Ai < A 2 < ■ ■ ■ < Afc+i; 

Let Ai = A + 7 / (as in 
Let Amin = niin{Ai, 7 }, and let r be its algebraic multiplicity; 
Dehne P^ ^ and c/y ^ P[g-, 

Compute ttj = {g\i)j for j = 1, k + 1 and 0^+2 = v^||fi'||i - 
if Amin > 0 and 0 ( 0 ) > 0 then 

cr* = 0 and compute p* from ( [IT| with r* = 7 ; 
else if Amin < 0 and 0(—Amin) > 0 then 

rr* — —\ 

^ — /'miri) _ 

Compute p* using ( [I^ ; 
if A min < 0 then 

Compute 2 ;* using (14); 
p* ^p* + z*] 


2 . 

2 ) 


end 

else 

Use Newton’s method to hnd a*, a root of 0, in (max{—Amin, 0}, 00 ); 
Compute p* from ( [II] ) with r* = a* + 7 ; 

end 


ALGORITHM 1: Orthonormal Basis SRI method 
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3.1. Newton’s method. Newton’s method is used to find a root of 0 ((t) whenever 

1 


lim 0 (cr) = lim 


1 


Since ||p(cr)|| does not exist at the eigenvalues of B, we first define the continuous 
extension of 0(a), whose domain is all of 3?. Let = {g\\)i for 1 < i < k + 1, 
ak+2 = llfi'alli and Xk+2 = 7 - Combining the terms in (10) that correspond to the 


same eigenvalues and eliminating all terms with zero numerators, we have that for 
lb('^)lP can be written as 


k +2 


IbWf = E 




E 


such that for z = 1, ...,£, Oj 7 ^ 0 and Aj are distinct eigenvalues of B with Ai < 
A 2 < ■ ■ ■ < A^. Note that the last sum is well-defined for a = Xj ^ Aj, for 1 < z < £. 
Then, the continuous extension 0(a) of 0(a) is given by: 


0(a) 


f 1 

“5 


E 


^717 + -)= 


if a = —Aj, 1 < z < £ 

-- otherwise. 

0 


A crucial characteristic of 0 is that it takes on the value of the limit of 0 at a = —A*, 
for 1 < z < A; -|- 2. In other words, for each z G {1 ,... ,k + 2}, 

lim 0(a) = 0(-Aj). 

cr—>■ —Ai 


The derivative of 0(a) is used only for Newton’s method and is computed as follows: 


a = 


E 


a? 


3 

'2 i 


E 


Cif 


^(Aj-Fa)2y (A* + a) 


if a 7 ^ —Aj, 1 < i < £. 


(17) 


Note that 0'(—Aj) exists as long as —Aj —Aj, for 1 < z < £. Furthermore, 
for a > —Ai, 0^(a) > 0, i.e., 0(a) is strictly increasing on the interval [—Ai,oo). 
Finally, it can be shown that (p'^cr) < 0 for a > —Ai, i.e., 0(a) is concave on the 
interval [—Ai,oo). For illustrative purposes, we plot examples of 0(a) in Fig. 1 
for the different cases we considered in this Section [H Note that we use Newton’s 
method to find a* when (a) Amin > 0 and 0(0) (see Figs. 1(b) and (c)), or (b) 
Amin < 0 and 0(—Amin) < 0 (see Figs. 1(d) and (e)). 

We now define an initial iterate such that Newton’s method is guaranteed to 
converge to a* monotonically. 


Theorem 3. Suppose 0(max{O, —Amin}) < 0. Let 

a = max | ^ — Ajl = — Aj (18) 

i<i<k +2 ( 5 J (5 

for some 1 < j < k + 2. Newton’s method applied to 0(a) with initial iterate 
0 -+) A max{ 0 , d} is guaranteed to converge to a* monotonically. 














10 


J. BRUST, J. B. ERWAY, AND R. F. MARCIA 



Figure 1. Graphs of the function 0(cr). (a) The positive-dehnite 
case where the unconstrained minimizer is within the trust-region 
radius, i.e., 0(0) > 0, and a* = 0. (b) The positive-dehnite case 
where the unconstrained minimizer is infeasible, i.e., 0(0) < 0. (c) 
The singular case where Ai = Amin = 0. (d) The indehnite case where 
Ai = Amin < 0. (e) When the coefficients a, corresponding to Amin are 
all 0, 0((j) does not have a singularity at Amin- Note that this case is 
not the hard case since 0(—Amin) < 0. (f) The hard case where there 
does not exist a* > —Amin such that (^(cr*) = 0. 
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Proof. Since 0(cr) is strictly increasing and concave on [—Amm, cxd) and 0(cr*) = 0, 
it is sufficient to show that (i) —Amin < < cr*, and (ii) exists (see e.g., 

ini)- 

We note that a > —Amin, and thus, > maxjO, —Amin} > —Amin- To show 
that (j® < cr*, we show that < 0 (cr*) = 0 . 

If d = \aj\/5 — Xj with \aj\ ^ 0, then evaluating ||p(cr)|| at a = d yields 


k+2 




> 


i=i (Aj + af 


(A, + M 

and thus, 0(d) < 0. Since 0(max{O, — Amin}) < 0, then 0(cr(°)) < 0. If \aj\ = 
0, then d = —Xj. Since -A* < —Amin for all i, a = —Amin- Thus, = 

0(niax{O, —Amin}) < 0. Consequently, 0(cr(°^) < 0, and therefore, < a* since 
0 ((t) is monotonically increasing. 

Next, we show that 0'(cr(°)) exists. On the interval (—Amin, oo), ^{o') is differen¬ 
tiable (see (17)). Therefore, if > —Amin, then 0'((T^°^) exists. If <7^°^ = —A min , 
, which implies that Oi = 


= <52 


then d = —A 
From the dehnition of 0(cr), A 
cr = -Amin = Ct(°A □ 


= Oj, = 0 or ak +2 = 0 (see (18)). 
7 ^ Aj for 1 < i < i. Thus, 0(cr) is differentiable at 


□ 


We note that when 7 ^ 0 in (18), d is the largest a that solves the secular 
equation with the inhnity norm: 


)(< 5 -) = 


n cr 




We illustrate the choice of initial iterate for Newton’s method in Fig. 2. 




(b) 


Figure 2. Choice of initial iterate for Newton’s method, (a) If aj 7 ^ 
0 in (18), then d corresponds to the largest root of 0 oo(cr) (in red). 
Here, —Amin > 0, and therefore cr(°) = d. (b) If aj = 0 in (18), then 
Amin 7 ^ Ai, and therefore, 0 (cr) is differentiable at —Amin since 0 (cr) is 
differentiable on (—Ai,cx)). Here, —Amin > 0, and thus, cr^^^ = d = 
Amin - 


Finally, we present Newton’s method for computing cr*. 
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Define tolerance r > 0; 
if 0(niax{O, —Amin}) < 0 then 
a = niaxi<j<fc+2 ^ - A^; 
a = maxjO, a}; 
while |0(cr)| > r do 
a = a - 0(a)/0'(cr); 
end 
a* = a; 

else if Amin < 0 then 

n-* — — X 
^ 

else 

a* = 0; 

end 

ALGORITHM 2: Newton’s method for compnting a* 


4. Numerical experiments 

In this section, we demonstrate the accnracy of the proposed OBS algorithm 
implemented in MATLAB to solve limited-memory SRI trust-region subproblems. 
For the experiments, hve sets of experiments composed of problems of various sizes 
were generated using random data. The Newton method to hnd a root of 0 was 
terminated when the ith iterate satished < ||0 (ct*'°^)|| -1- r, where 

denotes the initial iterate for Newton’s method and r = 1.0 x 10“^°. This is the 
only stopping criteria used by the OBS method since other aspects of this method 
compute solutions by formula. The problem sizes n range from n = 10^ to n = 10^. 
The number of limited-memory updates k was set to 4, and thus k + 1 = 5, and 
7 = 0.5 unless otherwise specihed below. The pairs S and F, both n x {k + 1) 
matrices, were generated from random data. Finally, g was generated by random 
data unless otherwise stated. The hve sets of experiments are intended to be 
comprehensive: They include the unconstrained case and the three cases discussed 
in Section The hve experiments are as follows: 

(1) The matrix B is positive dehnite with ||ptj|| < S: We ensure 4/ and M are 
such that B is strictly positive dehnite by altering the spectral decomposi¬ 
tion of RMR^. We choose S = /i||pn||, where fi = 1.25, to guarantee that 
the unconstrained minimizer is feasible. The graph of 0(cr) corresponding 
to this case is illustrated in Fig. 1(a). 

(2) The matrix B is positive dehnite with ||p„|| > S: We ensure 4/ and M are 
such that B is strictly positive dehnite by altering the spectral decompo¬ 
sition of RMR^. We choose 5 = /i||pn||, where /i is randomly generated 
between 0 and 1, to guarantee that the unconstrained minimizer is infea¬ 
sible. The graph of (^(cr) corresponding to this case is illustrated in Fig. 
1(b). 

(3) The matrix B is positive semidehnite and singular with p = —B^g infeasible: 
We ensure 4/ and M are such that B is positive semidehnite and singular 
by altering the spectral decomposition of RMR^. Two cases are tested: (a) 
0(0) < 0 and (b) 0(0) > 0. Case (a) occurs when 5 = (1 -|- /r)||p„||, where 
/i is randomly generated between 0 and 1; case (b) occurs when S = /i||Pm||, 
where p is randomly generated between 0 and 1. The graph of 0(cr) in case 
(a) corresponds to Fig. 1(c). In case (b), a* = 0 for i = 1,..., r, and thus. 
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0 (cr) does not have a singularity at a = 0 , implying the graph of 0 ((t) for 
this case corresponds to Fig i(a). 

(4) The matrix B is indehnite with 0(—Amin) < 0: We ensure T and M are such 

that B is indehnite by altering the spectral decomposition of RMR^. We 
test two subcases: (a) the vector g is generated randomly, and (b) a random 
vector g is projected onto the orthogonal complement of F||^ G so that 
Qi = 0,i = where r = 2. For case (b), 8 = /i||pu||, where /x is 

randomly generated between 0 and 1, so that 0(—Amin) < 0. The graph of 
0(cr) in case (a) corresponds to Fig. 1(d), and 0(cr) in case (b) corresponds 
to Fig. 1(e). 

(5) The hard case {B is indehnite): We ensure T and M are such that B is 
indehnite by altering the spectral decomposition of RMR^. We test two 
subcases: (a) Amin = Ai = Ai + 7 < 0, and (b) Amin = 7 < 0. In both 
cases, 5 = (1 + /i)||pn||, where fi is randomly generated between 0 and 1 , 
so that 0(—Amin) > 0. The graph of 0 ((t) for both cases of the hard case 
corresponds to Fig. 1(f). 

We report the following: (1) opt 1 (abs) = ||(i? + cT*/)p*+ 5 f||, which corresponds 
to the norm of the error in the hrst optimality conditions; ( 2 ) opt 1 (rel) =(||(i? + 
a*I)p* + 5 '||)/|| 5 '||, which corresponds to the norm of the relative error in the hrst 
optimality conditions; (3) opt 2=a*\p* — 5|, which corresponds to the absolute 
error in the second optimality conditions; (4) |(;/)(cr*)|, which measures how well 
the secular equation is satished; and (5) Time. We ran each experiment hve times 
and report one representative result for each experiment. We show in Fig. 3 the 
computational time for each of the hve runs in each experiment. 

For comparison, we report results for the OBS method as well as the LSTRS 
method [201 EH- The LSTRS method solves large trust-region subproblems by 
converting the subproblems into parametrized eigenvalue problems. This method 
uses only matrix-vector products. For these tests, we suppressed all run-time output 
of the LSTRS method and supplied a routine to compute matrix-vector products 
using the factors in the compact formulation (see ([^), i.e., given a vector v, the 
product with B is computed as Bv ■<— 7 n-|-T(M(T^n)). Note that the computations 
of M and T are not included in the time counts for LSTRS. 

Table 1. Experiment 1: OBS method with B is positive dehnite 
and IIp^II < 8. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

3.24e-15 

1.03e-16 

0.OOe+00 

0.OOe+00 

2 .12e-02 

l.Oe+04 

1.21e-14 

1.21e-16 

0.OOe+00 

0.OOe+00 

2.76e-02 

l.Oe+05 

4.61e-14 

1.46e-16 

0.OOe+00 

0.OOe+00 

5.46e-02 

l.Oe+06 

1.08e-13 

1.08e-16 

0.OOe+00 

0.OOe+00 

5.34e-01 

l.Oe+07 

5.31e-13 

1.68e-16 

0.OOe+00 

0.OOe+00 

5.34e+00 


Tables 1 and 2 shows the results of Experiment 1. In all cases, the OBS method 
and the LSTRS method found global solutions of the trust-region subproblems. The 
relative error in the OBS method is smaller than the relative error in the LSTRS 
method. Moreover, the OBS method solved each subproblem in less time than the 
LSTRS method. 
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Table 2. Experiment 1: LSTRS method with B is positive dehnite 
and ||p„|| < 5. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

2.lle-05 

6.70e-07 

0.OOe+00 

0.OOe+00 

4.72e-01 

l.Oe+04 

8.27e-07 

8.28e-09 

0.OOe+00 

0.OOe+00 

4.98e-01 

l.Oe+05 

2.64e-07 

8.37e-10 

0.OOe+00 

0.OOe+00 

9.15e-01 

l.Oe+06 

3.54e-09 

3.53e-12 

0.OOe+00 

0.OOe+00 

7.08e+00 

l.Oe+07 

2.79e-09 

8.81e-13 

0.OOe+00 

0.OOe+00 

6.66e+01 


Table 3. Experiment 2: OBS method with B is positive dehnite 
and IIp^iII > 5. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

3.44e-15 

1.06e-16 

1.75e-09 

4.82e+01 

2.83e-02 

l.Oe+04 

1.35e-14 

1.35e-16 

5.83e-13 

1.99e+01 

2.70e-02 

l.Oe+05 

3.34e-14 

1.06e-16 

6.15e-13 

1.57e+01 

6.39e-02 

l.Oe+06 

9.58e-14 

9.58e-17 

1.30e-ll 

7.06e+01 

5.38e-01 

l.Oe+07 

4.49e-13 

1.42e-16 

5.39e-06 

1.08e+00 

5.37e+00 


Table 4. Experiment 2: LSTRS method with B is positive dehnite 
and llptill > S. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

1.32e-14 

4.05e-16 

6.25e-04 

4.82e+01 

4.44e-01 

l.Oe+04 

1.20e-13 

1.20e-15 

1.20e-03 

1.99e+01 

4.80e-01 

l.Oe+05 

5.45e-ll 

1.73e-13 

4.90e-04 

1.57e+01 

7.30e-01 

l.Oe+06 

4.68e-10 

4.68e-13 

1.35e-06 

7.06e+01 

4.56e+00 

l.Oe+07 

4.15e-05 

1.31e-08 

4.47e-05 

1.08e+00 

4.21e+01 


Tables 3 and 4 show the results of Experiment 2. In this case, the unconstrained 
minimizer is not inside the trust region, making the value of a* strictly positive. 
As in the hrst experiment, the OBS method appears to obtain solutions to higher 
accuracy (columns 1,2, and 3) and in less time (column 4) than the LSTRS method. 
Finally, it is worth noting that as n increases, the accuracy of the solutions obtained 
by the LSTRS method appears to degrade. 

Table 5. Experiment 3(a): OBS method with B is positive semi- 
dehnite and singular with ||i?^5'|| > 5. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

2.80e-14 

8.89e-16 

6.25e-10 

3.38e-01 

2.70e-02 

l.Oe+04 

1.17e-13 

1.16e-15 

1.18e-08 

1.03e-01 

3.36e-02 

l.Oe+05 

3.48e-12 

l.lOe-14 

2.16e-07 

8.75e-03 

6.43e-02 

l.Oe+06 

1.44e-ll 

1.44e-14 

1.48e-09 

3.62e-03 

5.44e-01 

l.Oe+07 

5.52e-10 

1.74e-13 

8.96e-09 

2.88e-03 

5.39e+00 


Tables 5 and 6 display the results of Experiment 3(a). This is experiment is the 
hrst of two in which B is highly ill-conditioned. In this experiment, the LSTRS 
method appears unable to obtain solutions to high absolute accuracy (see column 
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Table 6. Experiment 3(a): LSTRS method with B is positive 
semidehnite and singular with ||-B^5'|| > 5. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

9.75e-03 

3.10e-04 

1.51e-16 

3.41e-01 

4.78e-01 

l.Oe+04 

7.93e-02 

7.91e-04 

2.65e-15 

1.07e-01 

5.69e-01 

l.Oe+05 

1.85e-01 

5.84e-04 

8.16e-16 

9.57e-03 

1.56e+00 

l.Oe+06 

1.29e-01 

1.29e-04 

6.04e-16 

1.70e-03 

1.28e+01 

l.Oe+07 

2.24e+03 

7.09e-01 

1.05e-10 

1.30e-06 

6.39e+01 


2 in Table 6). Moreover, the time required by the LSTRS to obtain solutions is, 
in some cases, signihcantly more than the time required by the OBS method. In 
contrast, the OBS method is able to obtain high accuracy solutions. Notice that 
the optimal values a* found by both methods appear to differ. Global solutions 
to the subproblems solved in Experiment 3(a) he on the boundary of the trust 
region. Because LSTRS was able to satisfy the second optimality condition to high 
accuracy but not the hrst, this suggests LSTRS’s solution p* lies on the boundary 
but there is some error in this solution. As n increases, the solution quality of the 
LSTRS method appears to decline with signihcant error in the case of n = 10^. 
In this experiment, the OBS method appears to hnd solutions to high accuracy in 
comparable time to other experiments; in contrast, the LSTRS method appears to 
have difficulty Ending global solutions. 

Table 7. Experiment 3(b); OBS method with B is positive semi- 

definite and singular with 115^(711 < S. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

4.lOe-15 

1.34e-16 

9.05e-10 

4.85e+01 

3.01e-02 

l.Oe+04 

1.Ole-14 

1.02e-16 

1.34e-ll 

6.98e+00 

4.36e-02 

l.Oe+05 

3.03e-14 

9.55e-17 

7.99e-14 

2.25e+01 

6.70e-02 

l.Oe+06 

1.39e-13 

1.39e-16 

4.18e-12 

3.42e+00 

5.41e-01 

l.Oe+07 

3.46e-13 

1.09e-16 

1.28e-ll 

1.08e+00 

5.37e+00 


Table 8. Experiment 3(b): LSTRS method with B is positive 
semidehnite and singular with ||i?^(7|| < 5. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

9.40e-15 

2.97e-16 

8.19e-04 

4.85e+01 

4.42e-01 

l.Oe+04 

2.06e-12 

2.07e-14 

6.59e-04 

6.98e+00 

4.79e-01 

l.Oe+05 

1.69e-ll 

5.34e-14 

4.27e-05 

2.25e+01 

7.43e-01 

l.Oe+06 

6.27e-08 

6.28e-ll 

6.19e-05 

3.42e+00 

4.60e+00 

l.Oe+07 

4.28e-05 

1.35e-08 

2.59e-05 

1.08e+00 

6.29e+01 


The results for Experiment 3(b) are shown in Tables 7 and 8. This is the second 
experiment involving ill-conditioned matrices. As with Experiment 3(a), the OBS 
method is able to obtain high-accuracy solutions in generally less time than the 
LSTRS method. The accuracy obtained by the LSTRS method appears to degrade 
as the size of the problem increases. In this experiment, the global solution always 
lies on the boundary, but the larger residuals associated the second optimality 
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condition in Table 8 indicate that the computed solutions by LSTRS do not lie on 
the boundary. 

Table 9. Experiment 4(a): OBS method with B is indehnite with 
Amin) < 0. The vector g is randomly generated. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

2.83e-15 

9.04e-17 

3.57e-12 

1.89e+02 

3.05e-02 

l.Oe+04 

1.27e-14 

1.27e-16 

1.53e-09 

1.18e+02 

3.99e-02 

l.Oe+05 

3.42e-14 

1.08e-16 

9.15e-13 

3.92e+02 

6.40e-02 

l.Oe+06 

1.19e-13 

1.20e-16 

4.79e-12 

5.39e+03 

5.43e-01 

l.Oe+07 

3.46e-13 

1.09e-16 

8.18e-ll 

1.94e+04 

5.35e+00 


Table 10. Experiment 4(a): LSTRS method with B is indehnite 
with 0(—Amin) < 0. The vector g is randomly generated. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

4.92e-14 

1.57e-15 

5.40e-04 

1.89e+02 

4.40e-01 

l.Oe+04 

2.82e-14 

2.79e-16 

1.03e-03 

1.18e+02 

4.80e-01 

l.Oe+05 

2.lle-13 

6.69e-16 

2.68e-06 

3.92e+02 

7.24e-01 

l.Oe+06 

2.93e-ll 

2.94e-14 

1.38e-07 

5.39e+03 

4.49e+00 

l.Oe+07 

1.81e-10 

5.74e-14 

3.19e-10 

1.94e+04 

4.12e+01 


The results for Experiment 4(a) are displayed in Tables 9 and 10. Both methods 
found solutions that satished the hrst optimality conditions to high accuracy. The 
overall solution quality from the OBS method appears better in the sense that the 
residuals for both optimality conditions in Table 9 are smaller than the residuals 
for both optimality conditions in Table 10. Finally, the OBS method took less time 
to solve the subproblem than the LSTRS method. 

Table 11. Experiment 4(b): OBS method with B is indehnite with 
0(—Amin) < 0. The vector g lies in the orthogonal complement of 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

3.42e-15 

1.07e-16 

1.17e-09 

1.31e+01 

2.91e-02 

l.Oe+04 

1.38e-14 

1.38e-16 

1.50e-14 

2.81e+00 

3.16e-02 

l.Oe+05 

3.17e-14 

l.OOe-16 

3.55e-13 

1.82e+01 

6.66e-02 

l.Oe+06 

1.30e-13 

1.30e-16 

1.76e-12 

4.76e+00 

5.46e-01 

l.Oe+07 

3.14e-13 

9.94e-17 

4.36e-ll 

7.58e+01 

5.36e+00 


The results of Experiment 4(b) are in Tables 11 and 12. Both methods solved 
the subproblem to high accuracy as measured by the hrst optimality condition; 
however, the OBS method solved the subproblem to signihcantly better accuracy as 
measured by the second optimality condition than the LSTRS method. All residual 
associated with the hrst and second optimality condition are less for the solution 
obtained by the OBS method. Moreover, the time required to hud solutions was 
less for the OBS method. 
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Table 12. Experiment 4(b): LSTRS method with B is indehnite 
with 0(—Amin) < 0. The vector g lies in the orthogonal complement 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

1.16e-14 

3.64e-16 

1.24e-03 

1.31e+01 

4.42e-01 

l.Oe+04 

2.48e-12 

2.49e-14 

1.02e-04 

2.81e+00 

4.70e-01 

l.Oe+05 

1.50e-10 

4.75e-13 

2.82e-04 

1.82e+01 

7.30e-01 

l.Oe+06 

1.65e-08 

1.65e-ll 

9.70e-05 

4.76e+00 

4.65e+00 

l.Oe+07 

2.08e-07 

6.58e-ll 

1.06e-05 

7.58e+01 

4.21e+01 


Table 13. Experiment 5(a): The OBS method in the hard case {B 
is indehnite) and Amin = Ai = Ai + 7 < 0. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

1.29e-14 

4.34e-16 

1.93e-16 

4.35e-01 

3.38e-02 

l.Oe+04 

5.87e-14 

5.86e-16 

2.59e-14 

6.08e-01 

2.73e-02 

l.Oe+05 

2.34e-12 

7.43e-15 

5.79e-14 

8.15e+00 

8.08e-02 

l.Oe+06 

1.33e-ll 

1.33e-14 

1.19e-12 

3.97e+00 

6.72e-01 

l.Oe+07 

1.67e-10 

5.28e-14 

4.43e-12 

5.27e-01 

6.71e+00 


Table 14. Experiment 5(a): The LSTRS method in the hard case 
{B is indehnite) and Amin = Ai = Ai + 7 < 0. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

2.lOe-05 

7.07e-07 

1.16e-15 

4.35e-01 

4.70e-01 

l.Oe+04 

3.88e+00 

3.87e-02 

1.50e-03 

6.08e-01 

4.71e-01 

l.Oe+05 

1.27e+02 

4.01e-01 

5.72e-04 

8.15e+00 

7.65e-01 

l.Oe+06 

2.04e+02 

2.04e-01 

1.45e-04 

3.97e+00 

4.59e+00 

l.Oe+07 

1.64e+03 

5.17e-01 

2.30e-05 

5.27e-01 

4.23e+01 


In the hard case with Amin being a nontrivial eigenvalue, the OBS method was 
obtain global solutions to the subproblems; however, the LSTRS had difficulty hnd- 
ing high-accuracy solutions for all problem sizes. In particular, as n increases, the 
solution quality of the LSTRS method appears to decline with signihcant error in 
the case of n = 10^. In all cases, the time required by the OBS method to hnd a 
solution was less than that of the time required by the LSTRS method. 

Table 15. Experiment 5(b): The OBS method in the hard case {B 
is indehnite) and Amin = 7 < 0. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

3.52e-15 

l.lle-16 

3.53e-09 

6.35e+01 

2.93e-02 

l.Oe+04 

9.50e-15 

9.48e-17 

1.16e-14 

2.10e+02 

3.82e-02 

l.Oe+05 

3.Ole-14 

9.50e-17 

4.49e-13 

4.49e+02 

6.71e-02 

l.Oe+06 

9.48e-14 

9.47e-17 

6.86e-12 

1.34e+04 

5.32e-01 

l.Oe+07 

3.40e-13 

1.07e-16 

2.97e-12 

8.91e+03 

5.36e+00 
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Table 16. Experiment 5(b): The LSTRS method in the hard case 
{B is indehnite) and Amin = 7 < 0. 


n 

opt 1 (abs) 

opt 1 (rel) 

opt 2 

a* 

Time 

l.Oe+03 

2.24e-14 

7.12e-16 

7.36e-04 

6.35e+01 

4.41e-01 

l.Oe+04 

6.35e-14 

6.33e-16 

1.92e-06 

2.lOe+02 

5.02e-01 

l.Oe+05 

2.26e-13 

7.14e-16 

5.09e-08 

4.49e+02 

7.49e-01 

l.Oe+06 

6.61e-12 

6.61e-15 

4.76e-08 

1.34e+04 

4.32e+00 

l.Oe+07 

8.77e-ll 

2.77e-14 

1.05e-08 

8.91e+03 

4.09e+01 


The results of Experiment 5(b) are in Tables 15 and 16. Unlike in Experiment 
5(a), the LSTRS method was able to hnd solutions to high accuracy. In all cases, 
the OBS method was able to hnd solutions with higher accuracy than the LSTRS 
method and in less time. 


5. Concluding remarks 

In this paper, we presented the OBS method, which solves trust-region subprob¬ 
lems of the form ([^ where B is a large L-SRl matrix. The OBS method uses two 
main strategies. In one strategy, a* is computed from Newton’s method and initial¬ 
ized at a point where Newton’s method is guaranteed to converge monotonically to 
a*. With a* in hand, p* is computed directly by formula. For the other strategy, we 
propose a method that relies on an orthonormal basis to directly compute p*. (In 
this case, a* can be determined from the spectral decomposition of B.) Numerical 
experiments suggest that the OBS method is able to solve large L-SRl trust-region 
subproblems to high accuracy. Moreover, the method appears to be more robust 
than the LSTRS method, which does not exploit the specihc structure of B. In 
particular, the proposed OBS method achieves high accuracy in less time in all of 
the experiments and in all measures of optimality than the LSTRS method. Future 
research will consider the best implementation of the OBS method in a trust-region 
method (see, for example, m), including initialization of 7 and rules for updating 
the matrices S and Y containing the quasi-Newton pairs. 
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Figure 3. Semi-log plots of the computational times (in seconds). 
Each experiment was run Eve times; computational time for the 
LSTRS and OBS method are shown for each run. In all cases, the 
OBS method outperforms LSTRS in terms of computational time. 
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